Power Spectra of the Total Occupancy in the Totally Asymmetric Simple Exclusion 
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As a solvable and broadly applicable model system, the totally asymmetric exclusion process 
enjoys iconic status in the theory of non-equilibrium phase transitions. Here, we focus on the 
time dependence of the total number of particles on a 1-dimensional open lattice, and its power 
spectrum. Using both Monte Carlo simulations and analytic methods, we explore its behavior in 
different characteristic regimes. In the maximal current phase and on the coexistence line (between 
high/low density phases), the power spectrum displays algebraic decay, with exponents —1.62 and 
—2.00, respectively. Deep within the high/low density phases, we find pronounced oscillations, 
which damp into power laws. This behavior can be understood in terms of driven biased diffusion 
with conserved noise in the bulk. 

PACS numbers: 02.50.Ey, 05.60.-k, 05.40.-a, 87.10.+C 
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Introduction. The collective behavior of many con- 
stituents, interacting under non- equilibrium conditions, 
is far from well understood. Yet, such systems are ubiq- 
uitous in nature, from molecular biology at the nanoscale 
to infractructure networks at the global level. In physics, 
attacks on such highly complex systems often begin with 
seemingly small steps, defining simple models that are 
both tractable and believed to capture the essentials of 
the original problem. On this stage, the totally asymmet- 
ric simple exclusion process (TASEP) [1] plays a key role, 
much like the Ising model in equilibrium statistical me- 
chanics. With open boundaries, particles stochastically 
enter a one-lane lattice from one end, hop to the next 
site if it is empty, and leave at the opposite end. This 
simple transport model provides the first crucial steps 
towards the modeling of realistic processes, such as pro- 
tein synthesis, surface growth, and vehicular traffic Q. 
Further, it generates much theoretical interest: Despite 
its simplicity, its stationary state displays several phases 
and interesting algebraic correlations, much of which is 
known analytically 

In this study, we focus on another simple quantity in 
TASEP: the total number of particles in the system at 
time t, N(i). In the steady state, its time average should 
be the ensemble average, which can be easily computed 
from the known stationary distribution. However, as a 
fluctuating quantity, its power spectrum, I{uS), contains 
time-correlation information which is not as easily acces- 
sible. Although such correlations have been investigated 
recently [3] , we find interesting behavior undetected pre- 
viously: oscillatory behavior in I(oS) for the high and 
low density phases. While the previous study concerns 
correlations within the bulk of an infinite system, the 
new feature here is that our / (ui) carries information on 
the entirety of a finite system. Since many physical sys- 
tems are far from the thermodynamic limit (e.g., mRNAs 
containing around 10 3 codons or fewer), such finite-size 
effects can be physically significant. In the remainder of 



this letter, we briefly summarize the details of the model, 
our simulation methods, report our findings, and provide 
theoretical explanations for the phenomena. 

The model and simulation results. A standard TASEP 
consists of particles hopping on a lattice of length L with 
site label i. If the first site is empty, a particle enters 
the system with rate a. A particle on the last site leaves 
the system with rate [3. In the bulk, a particle always 
hops onto the next site (with unit rate) if the site is 
empty; otherwise, it remains stationary. A configuration 
at any particular time is described by a set of occupation 
numbers {nt} (n = 0, 1 for a vacant vs. occupied site). 
In all our simulations, we start with an empty lattice and 
use random sequential updates, i.e., in each Monte Carlo 
step (MCS), we make N + 1 random attempts to move 
a particle. Since the system is stochastic, a complete 
description requires P({ni},t), the probability to find 
the system in configuration {n^} at time t. The evolution 
of P is governed by a simple master equation. Though 
linear and typically easy to write, this equation cannot 
be solved in general. Yet, in the long time limit, P settles 
into a t-independent distribution: P*, the exact form of 
which is known [3|. 

The macroscopic properties of these stationary states 
can be categorized in terms of three distinct phases, con- 
veniently displayed as a phase diagram (Fig. QJi) in the 
a- (3 plane. Thanks to the exact solution, P* , many prop- 
erties can be computed analytically. The three regions 
are associated with the high density (HD), low density 
(LD), and maximum current (MC) phases. They are dis- 
tinguished by their average local densities, pi = (rii/L), 
where (•) is an average over the stationary P* ({n^}). 
In the thermodynamic limit and deep in the bulk, pi is 
p + = 1 — (3, p- = a, and 1/2 for, respectively, the HD, 
LD, and MC phases. All aspects of the LD and HD 
phases can be directly related through particle-hole sym- 
metry, so that we will concentrate only on the LD phase 
in this letter. For the LD phase, the deviations from p_ 
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FIG. 1: (a) Phase diagram of TASEP. On the dashed line, 
high and low density domains coexist, with a microscopic in- 
terface, known as the "shock." (b) Typical time series (up 
to 600K MCS shown) of N(t), for L = 1000, in the three 
distinct regimes: coexistence, MC and LD (color online). 

are confined to a microscopic boundary layer near the 
exit. In the MC phase, the profile decays algebraically 
into the bulk, while the system displays behavior com- 
parable to critical phenomena in equilibrium statistical 
mechanics. The transitions at the boundaries of the MC 
phase are continuous. By contrast, a discontinuity oc- 
curs across the HD/LD line. On this line (a = (3 < 1/2), 
typical configurations display the coexistence of regions 
with p + and p_, joined by a microscopically sharp in- 
terface, known as a "shock". The shock is delocalized, 
performing a random walk over the entire lattice, so that 
the average pi is linear over the interval [a, 1 — a]. 

Our main interest here is N(t) = anc ^ ^s as- 

sociated power spectrum: I(lu). Though its time average 
(in the steady state) is just (rii) = pL, I(lu) cannot 
be accessed from P*. Specifically, we take a measure- 
ment of N every 100 MCS (after discarding the first 100K 
MCS to allow the system to reach steady state) and la- 
bel these by t = 1,2, ...,T. Typical N(t)'s in the three 
regimes are shown in Fig. Q}). For all but a few of our 
runs, T = 10 4 (10 6 MCS), so that we can define a Fourier 
transform: N (uj) = J2t=i N(t)e luJt , where uj = 2-nrajT 
with m = 0, 1, ... . To obtain the average power spec- 
trum, we carry out typically 100 such runs and construct 



N(uj) 



(1) 



The important control parameters for this system are L, 
a, and (3. We investigated a range of L's from 250 to 
32000. Here, we will show the results for L = 1000, as 
well as some for the largest size. Also, we will present 
mainly data for three representative points in the phase 
diagram. Labeled by (a, ft), they are (0.7,0.7), (0.3,0.7), 
and (0.3,0.3), corresponding to, respectively, the MC, 
LD phases and the coexistence line. Systematic studies 
of the rest of the a-/3 square will be reported elsewhere 
1- 

In the inset of Fig. [2j we show the power spectra for 
the MC phase. Away from high and low uj values, / (uj) 
appears to obey a power law: / oc uj~ p . The saturation 
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FIG. 2: Power spectrum in the MC phase, a — /3 = 0.7. with 
L — 1000. The inset shows a log-log plot of I (u>). The main 
figure shows C = In (uj p I) vs. m = ujT '/2tt (color online). 

at small uj is due to finite L and the crossover scales with 
the expected L~ 3 / 2 (see below). At large uj, we are pre- 
sumably probing an ultraviolet cutoff, due to a discrete 
lattice and MCS. To expose the power p more effectively, 
we display uj p I in a log-log plot (Fig. [2J. The best fit 
seems to be p — 1.62 (line, blue online), though the data 
are still consistent with 5/3 (dots, red online), a value fa- 
vored by the theoretical considerations below. Whether 
the difference of ~ 3% is truly significant remains to be 
explored [B| . We have also simulated other points in the 
MC phase and find all data to be statistically indistin- 
guishable from the set here. The results for the coex- 
istence line, at a = /3 = 0.3, are very similar, except 
that they confirm the expected random walk behavior, 
with a power of 2. The most remarkable phenomena are 
found in the HD/LD phases. As illustrated in Fig. [3l 
I (uj) displays large, damped oscillations. The examples 
in the inset are especially striking: An L = 32K lattice is 
used here, with three a's. In the next section, we briefly 
present the theory which accounts for these properties. 

Theoretical approaches. To understand the behavior 
of the power spectra above, we rely on different (though 
related) approaches for the three regimes. The simplest 
case is the coexistence line, where the shock performs a 
random walk. As a result, N (t) is also just a random 
walk (Fig. [TJd) in the interval [p~L, p+L] with reflecting 
boundary conditions [ft]. The power spectrum of this 
process is well understood, displaying the uj~ 2 behavior. 
Of course, for small uj (and long times), / (uj) saturates 
to a finite value, controlled by the allowed interveral for 
N. Later we will show that the expected crossover ex- 
hibits data collapse with a scaling variable ujL 2 . The next 
regime, considered the most challenging theoretically, is 
the MC phase. Here, the associated non-equilibrium 
steady state is "critical", displaying (power law) corre- 
lations and anomalous exponents. At the coarse-grained 
level, our system is described by a continuous density, 
p(x,t), and may be studied as a stochastic field theory. 
Indeed, the TASEP is a driven diffusive system (jj in 
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FIG. 3: Power spectrum I (oj) vs. m = u)T/2n in the LD 
phase (a = 1 — j3 = 0.3) with L = 1000. The inset shows a 
similar plot for a, L — 32000 lattice, with a = 1— j3 = 0.2, 0.3, 
and 0.4. The units for log (I) are arbitrary, and the curves 
are displaced for clarity (color online). 

one dimension. Deviations of p about its average here, 
4>(x,t) = p(x,t) — 1/2, satisfy the noisy Burgers equa- 
tion [g| : d t (f> = i;dl4> + d x (f> 2 + d x i], where rj is a Gaussian 
noise. Within the bulk of an infinite system, the dynamic 
two-point function (<p (x, t) <p (x' , t')) is translationally in- 
variant and assumes a scaling form 0, 0, 0] : 



(2) 



where £ = x' — x and r = t' — t. If we naively con- 
sider J d£C (£, r) as a candidate for our power spectrum, 
we would find that Eqn. @ leads to J dsf (s), indepen- 
dent of t! This surprising conlusion is consistent with 
simply setting k — in the results of [13] . On closer 
examination, however, we should seek the Fourier trans- 
form of Jq dxdx' ((f) (x, 0) <p (x', t)) in a finite system. If 
we approximate this two-point function by C (£, t), we 
are led to d£ [2 (L — £) C (£, t)], which is of the form 
t 2 / 3 g (L/t 2 ^ 3 ) . Thus, its Fourier transform would take 
the scaling form / (oj) = oj~ 5 ' 3 G (Loj 2 ^) . Assuming that 
G approaches a positive constant for oj ^s> L~ 3 / 2 , this the- 
ory predicts I (oj) cx w~ 5 / 3 . As discussed above, our data 
are consistent with this power, although 1.62 appears to 
be a better fit. Though a 3% disagreement seems minor, 
we found similar discrepancies in other observables such 
as (N (t) N (i')). These issues, along with finite-size cor- 
rections to scaling, discrete space and time, and possible 
systematic errors, will be explored elsewhere [3]. 

Typically, the behavior deep in the HD and LD phases 
seems least interesting, with ordinary Gaussian fluctua- 
tions. Yet, / (oj) for this regime exhibits the most struc- 
ture. We now turn to a brief account of a simple, linear 
theory for the fluctuations that predicts these interesting 
oscillations, focusing on the LD phase only. 

Following standard routes for a continuum description 
[1,0] of a driven density p(x,t), we study small fluctua- 
tions, 4>, around p using the stochastic equation of mo- 



tion: d t (j)(x, t) = Dd 2 4> — vd x (f) + d x -q. Here, D is the (ef- 
fective) diffusion constant, v is the bias, and r\ is a Gaus- 
sian correlated (local) noise with zero mean. If we were 
to start from the microscopic {n{\ and take a naive con- 
tinuum limit, we would arrive at D = 1/2 and v = 1 — 2p 
(in units of lattice spacing and MCS). But we leave these 
as 0(1) parameters for now. The solution in Fourier 
space is easy: 4>{k, oj) = ikfj (k,oj) / \jJk 2 + ink — ioS\, 
where <j>{x,t) = J k u $ (k, oj) e^-"*) and J k = J dk/2ir. 
With N{t) = J^dxp{x,t) = pL + J^4>{x,t), we have 
N (oj) = f k (e ikL - 1) f\l [Dk 2 + ivk - ioj] for oj > 0. 
Writing (7~y (fc, oj) fj* (k\ oj 1 )) = AS (k - k') S (oj - oj') and 
replacing 5 (oj — oj) by T, we find the power spectrum to 
be 4^ J k | (e ikL - 1) / [Dk 2 + ivk - ioj] | 2 . Evaluting the 
integral, we arrive at 



I (oj) = = — Re 
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1 e lk - L - 1 



R(l-R) 2 R*(l + R*f 



(3) 



where R = y/l - AiDoj/v 2 and fc± = iv (-1 ± R) /2D. 

To shed light on this complex expression, we consider 
three frequency regimes: (a) oj — > 0, e.g., m = 1,2,3..., 
(b) intermediate oj's (oj <C v 2 /D, i.e., 1 <C m <C T), 
and (c) oj — > oo. In regime (a), I approaches a constant, 
which can be used to fit A. A further simplication is that 

the first term in Re [...] is O (j(vL/D) 2 ^j , while the second 

term is O (1). Thus, the latter can be largely ignored, if 
we take, e.g., L — 32K, as in the inset of Fig. [3l In 
regime (b), we keep O (oj 2 ) corrections in R and k±, i.e., 

2iDoj 2D 2 oj 2 oj iDoj 2 

1 7. 1 -. — and 

v v 

respectively. The result shows damped oscillations man- 
ifestly: 
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The minima are near integer multiples of 2irv/L in oj 
(i.e., vT/L in m). Remarkably, the minima in Fig. [3] 
are very close to these values if we naively substitute 
(0.4, 10 6 , 10 3 ) for (v, T, L). Encouraged by this approach, 
we attempt to fit the L = 32K data with the full Eqn. 
Q. The result for a = 0.3, shown in Fig. [4J is surpringly 
good, indicating that this simple theory has captured the 
essentials of our system. 

The physical origins of the oscillations are noteworthy. 
They can be traced to the time it takes a fluctuation to 
traverse the entire lattice. Indeed, if there is no exclusion 
and each particle travels ballistically at velocity u, then 
N (t) will consist of a random series of unit increments 
(due to a), correlated with unit drops at times L/u later 
(assuming (3 = 1). The associated I will be proportional 
to z~ 2 sin 2 z, with z = ojL/2u, with zeros at multiples of 
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FIG. 4: Power spectrum log (7) vs. 
with L = 32000 (color online). 



m for a = 1 — (3 = 0.3 



2-ku/L. The effect of diffusion, embodied in D, is just to 
smear out the minima and fill in the zeros. Returning to 
the fit above, we were surprised that it requires D ~ 20 
and A ~ 1/500, values that seem remarkably far from 
the naive O (1) levels. Work is in progress to understand 
why these parameters are so seriously "renormalized." 

Beyond the oscillations, there is another subtle fea- 
ture in this intermediate to regime. As Eqn. ([4]) shows, 
the oscillations are completely damped when to 2 exceeds 
v^/DL, and I settles into uj~ 2 . Though not shown, this 
power is confirmed in the L — 32K case. However, on 
further scrutiny, we see that R — » i/tJ in the truly asymp- 
totic regime (c): w ^> v 2 /AD. Then, Eqn. ([3]) predicts 
a uj~ 3 / 2 decay [41]! If this regime is reached before the 
oscillations are fully damped, then the lo" 2 behavior of 
the "intermediate regime" will lie hidden. This is indeed 
true for the L = IK case, where the oscillations die out 
at much larger m's and only the 3/2 power is observed. 

Outlook. In addition to a more detailed and system- 
atic Monte Carlo study (e.g., larger systems, longer runs, 
finite-size scaling, etc.), there are a number of issues wor- 
thy of further pursuit. The continuum theory represents 
a rough first step and should be refined to a version with 
discrete space/time. Most intimately connected to this 
work is a better understanding of the origin of the large 
"renormalization effects" on the diffusion coefficient, D. 
Beyond interests specific to the simple TASEP, we are 
motivated by its applicability to protein synthesis, where 
the lattice (particles) model the mRNA (ribosomes) 
To capture more of the process in vivo, we should in- 
clude the effects of (i) particles that occupy I > 1 sites, 
(ii) inhomogeneous hopping rates, (iii) finite reservoir of 
particles which can enter the lattice, and (iv) competi- 
tion with other TASEPs (mRNAs) for a finite pool of 
particles. Some of these issues have been considered 
though none focused on the total occupancy, a quan- 
tity surely of interest in the context of biology. Needless 
to say, many other generalizations - all well motivated 
by physical systems - come to mind, e.g., many parti- 
cle species, multiple lanes, and high dimensions. In the 



past, power spectra have served us well, providing valu- 
able perspectives into stochastic systems in general [12 1 
and driven system in particular [13[. We hope this work 
will spark further interest to study power spectra in novel 
non-equilibrium systems far beyond physics, such as bi- 
ology, social networks, and finance. 
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